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CO oxidation on Pd(100) at technologically relevant pressure conditions: 
A first-principles kinetic Monte Carlo study 
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The possible importance of oxide formation for the catalytic activity of transition metals in 
heterogenous oxidation catalysis has evoked a lively discussion over the recent years. On the more 
noble transition metals (like Pd, Pt or Ag) the low stability of the common bulk oxides suggests 
primarily sub-nanometer thin oxide films, so-called surface oxides, as potential candidates that 
may be stabilized under gas phase conditions representative of technological oxidation catalysis. 
We address this issue for the Pd(100) model catalyst surface with first-principles kinetic Monte 
Carlo (kMC) simulations that assess the stability of the well-characterized (a/5 x \/E)R27° surface 
oxide during steady-state CO oxidation. Our results show that at ambient pressure conditions 
the surface oxide is stabilized at the surface up to CO:02 partial pressure ratios just around the 
catalytically most relevant stoichiometric feeds (pco : Po 2 =2:1). The precise value depends 
sensitively on temperature, so that both local pressure and temperature fluctuations may induce 
a continuous formation and decomposition of oxidic phases during steady-state operation under 
ambient stoichiometric conditions. 

PACS numbers: 82.65.+r, 68.43.Bc, 68.43.De, 68.35.Md 
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I. INTRODUCTION 

The large difference in catalytic activity of "ruthe- 
nium" catalysts towards CO oxidation under ultra-high 
vacuum and ambient pressure conditions^ has recently 
been controversially discussed in terms of a possible ox- 
ide formation at the surface in technologically relevant 
environments^. The well known much lower thermody- 
namic stability of the bulk oxides of Pd or Pt£ suggests at 
first glance that this discussion has no direct relevance for 
these more noble metals that are equally, if not more of- 
ten used in oxidation catalysis. Using a constrained first- 
principles thermodynamics approach 6 we could show in 
a preceding study^, however, that the situation is not 
that clear-cut when considering for the possibility of 
sub- nanometer thin oxide films, so-called surface oxides. 
Mapping out a wide range of (T,po2,Pco)- conc htions 
we specifically found for the Pd(100) surface that the 
stability range of the well-characterized (a/5 x yfb)R27° 
(a/5 in brevity) surface oxide struct ure& & 1 1 1 1 1 1 2 does 
indeed extend up to pressure and temperature condi- 
tions relevant for technological CO oxidation catalysis 
(T = 300-600K, Po . 2 =p C o ~ latm). 

This insight motivates to center a detailed kinetic 
modeling of the surface structure and composition dur- 
ing steady-state CO oxidation at technologically relevant 
pressure conditions on the a/5 surface oxide. More pre- 
cisely, the question of a possible oxide formation at the 
surface can be narrowed down to addressing the stabil- 
ity, formation and decomposition of the a/5 surface oxide 
structure in the reactive environment. Theoretically, one 
is then still faced with the challenge to evaluate the kinet- 
ics of the surface elementary processes over time scales of 
the order of seconds or longer, while accurately account- 
ing for the evolving atomistic structure of the surface. 



In this study, we meet the time-scale challenge with a 
first-principles statistical mechanics approach, combin- 
ing density-functional theory (DFT) for an accurate de- 
scription of the kinetic parameters with kinetic Monte 
Carlo (kMC) simulations for the evaluation of the non- 
equilibrium statistical mechanics problem. At present, 
such an approach is preferentially carried out on the ba- 
sis of a lattice model for the surface in order to keep 
the number of required DFT-based parameters tractable. 
Due to the different periodicity and Pd-density of the a/5 
over layer ^ i 13 i 14 and the Pd(100) substrate, a combined 
lattice model embracing both structures will be highly 
involved and will require numerous kinetic parameters 
for a manifold of conceivable atomistic pathways for the 
complex structural transition between the oxidized and 
pristine Pd surface. Before embarking on such an explicit 
modeling of the formation and decomposition of the sur- 
face oxide structure in the simulations it is thus appro- 
priate to first assess whether the stability range of the 
v5 surface oxide does indeed extend to the catalytically 
relevant environments, when the kinetics of the ongoing 
CO2 formation is fully accounted for. In this study we 
determine this stability range by focusing on the onset 
of the decomposition of the surface oxide structure when 
the partial pressure ratio pco : Po 2 exceeds a critical 
value. We identify this onset once the average O concen- 
tration in the a/5 structure is reduced below the value 
corresponding to the intact surface oxide, which allows 
us to restrict our kMC modeling to the a/5 lattice. 

The central result of our kMC simulations, which we 
had briefly highlighted before^, is that at ambient pres- 
sures the surface oxide is stabilized at the surface up to 
C0:02 partial pressure ratios just around the catalyti- 
cally most relevant stoichiometric feeds (pco ■ Po 2 — 2 : 
1). This suggests an interpretation of the morphologi- 
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cal changes observed in recent reactor scanning-tunneling 
microscopy (STM) measurements 1 ^ under O-rich feeds in 
terms of a formation of the y/E structure at the surface. 
The precise value for the critical partial pressure ratio 
we compute depends furthermore sensitively on temper- 
ature, so that both local pressure and temperature fluctu- 
ations may induce a continuous formation and decompo- 
sition of the oxidic phase during steady-state operation 
under ambient stoichiometric conditions. We therefore 
expect that a full understanding of the catalytic activity 
of the Pd(100) surface in such technologically relevant 
environments can only be obtained from a modeling that 
explicitly accounts for the corresponding on-going struc- 
tural changes, which will then also be able to address 
the actual role played by the surface oxide for the cat- 
alytic activity itself. Nevertheless the rather high peak 
turnover frequencies obtained in our simulations even at 
the intact surface oxide show already that the formation 
of this structure in the reactive environment can not sim- 
plistically be equated with a "deactivation" of the model 
catalyst. 



II. THEORY 
A. First-principles kinetic Monte Carlo simulations 

Our target is to describe the steady-state catalytic 
CO oxidation over the \fb surface oxide by explicitly ac- 
counting for the kinetics of the different elementary pro- 
cesses taking place in the system. Properly evaluating 
the time evolution of such a system requires simulation 
cells that are large enough to capture the effects of cor- 
relation and spatial distribution of the chemicals at the 
surface. Furthermore, most of the considered processes 
are highly activated and occur on time scales that are 
orders of magnitude longer than e.g. a typical vibration 
(~ 10~ 12 s). Due to these so-called rare events we need 
to evaluate the statistical interplay between the different 
elementary processes over an extended time scale that 
can reach up to seconds or more. Since this is unfeasible 
for molecular dynamics simulations, we employ kinetic 
Monte Carlo simulations, where the time evolution is ef- 
ficiently coarse-grained to the rare-event dynamics. In a 
kMC simulation, each kMC-step corresponds to the ex- 
ecution of a rare event, while appropriately accounting 
for the short-time dynamics of the system in between. 
In the generated sequence of system states, each state 
is assumed to be independent of all previous states, i.e. 
the kMC algorithm provides an efficient numerical solu- 
tion to the Master equation of a continuous time Markov 
processiii^^^i. 

In order to keep the number of considered elementary 
processes in a tractable range the kMC simulations are 
performed on a lattice. Obviously, the chosen lattice 
model has to properly reflect the important structural 
features of the real system. For this model a list of all 
relevant elementary processes p on the lattice is set up 



with the corresponding rate constants k p that charac- 
terize the probability to escape from the present system 
state by the process p. Starting in a given initial system 
configuration the sum over the rate constants of all pos- 
sible processes in this current configuration is evaluated, 
fctot = Y] p kp- One process i is then randomly chosen by 
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k p < pifctot < 



p=0 



(1) 



with p\ g]0, 1] being a random number. Since the prob- 
ability of selecting a process i is weighted by its rate 
constant fcj a process with a large rate constant is more 
likely to be chosen than a process with a small rate con- 
stant. The selected process is executed and the system 
configuration is updated. This way the kMC algorithm 
simulates a sequence of Poisson processes and an explicit 
relationship between kMC time and real time can be es- 
tablished 21 . This evolution of real time after each kMC 
step is given by 



t->t- 



ln(p 2 ) 



(2) 



where p 2 €]0, 1] is a second random number—. 

The crucial ingredients to a meaningful kMC sim- 
ulation are therefore the identification of all relevant 
elementary processes p and the determination of the 
corresponding rate constants k p . For the latter we 
rely on the modern first-principles kinetic Monte Carlo 
approac h 22 ' 23 i 24 ' 25 i 26 ' 27 i 28 i 29 and calculate the rate con- 
stants using DFT and transition-state theory (TST). 
Concerning the elementary processes during the CO oxi- 
dation over the \/5 surface oxide we consider adsorption, 
desorption, diffusion and reaction of the chemicals, and 
follow the recipe of how to determine the rate constants 
for these processes from DFT and TST given in Ref. |29| . 
In the next Section the resulting expressions for the rate 
constants are briefly reviewed. The first-principles data 
required to evaluate the rate constants with these expres- 
sions is then given in Section [II El 



B. The rate constants 

1. Adsorption 

The rate constant of adsorption of a gas phase species 
i on a surface site of type st depends on the impingement 
rate and the local sticking coefficient, Si 
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Si, st (T) 



PiA 



V2 



7rmi«B 



fc R T 



(3) 



The only system-specific parameters entering the im- 
pingement rate are the mass m t of the impinging gas 
phase molecules and the surface area A. An explicit, ac- 
curate evaluation of the local sticking coefficient Si lS t is 
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highly involved and requires a determination of the full, 
high-dimensional potential energy surface (PES) of the 
adsorption process, on which molecular dynamics (MD) 
simulations have to be performed for a statistically rele- 
vant number of trajectories 3 ^. In the present study the 
local sticking coefficient is instead roughly estimated as 29 



A, 



A 



exp 



fc B r 



(4) 



diffusion pathway can be identified, harmonic TST can 
be applied to obtain the corresponding rate constant of 
diffusion 2 ^ 



fc^(T) = /^(T)(^) 



exp 



with 



rdiff.TST/rrrt 
J i,st— >st' V / 



AE^ 



i,st — *st' 



5 i,st — »at',TS 



(6) 
(7) 



where AEff t is the highest barrier along the minimum 
energy pathway (MEP) and Ai >st is the so-called active 
area. Within the concept of an active area it is assumed 
that only particles of a species i with an initial lateral po- 
sition within a certain area Aj iS t around the adsorption 
site st can actually stick to this site. The rate constant 
of adsorption in Eq. ([3]) will then actually become in- 
dependent of the total impingement area A and only the 
active area A iySt has to be known. f£*^ is a factor that re- 
duces the number of impinging gas phase particles by the 
fraction that is not traveling along the MEP and might 
thus be reflected at some higher barrier along a different 
pathway. 



and 



AEf 



El 



TTltOt 

&i,st 



(8) 



The diffusion barrier AE?f t ^ st> denotes the maximum 
barrier along the MEP of the diffusion process and is 
given by the energy difference of the transition state 

m%^ st >,Ts) and the initial state ( E l%) at site st - The 
reverse process of the diffusion st — > st' is simply the 
backward diffusion st' — > st and therefore the rate con- 
stants of these two processes have to fulfill the detailed 
balance criterion. To calculate the diffusion rate con- 

cliff T^ST 

stants the factor f i st l^ st , (given by the ratio of the vi- 
brational partition functions in the transition state (TS), 
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Zi7t-^st' TS' an d tn e initial state zj%) as well as the dif- 
fusion barriers have to be known. 



2. Desorption 



Since the desorption process is the time-reversed pro- 
cess of adsorption the desorption rate constants are con- 
nected to the adsorption rate constants by the detailed 
balance criterion and can thus be obtained by applying 29 
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TT'bind 



where AGi, s t(T,pi) is the change in the Gibbs free en- 
ergy between a particle in the gas phase and in the ad- 
sorbed state. AGi )S t(T,pi) is approximated by the dif- 
ference between the chemical potential of the particle 
in the gas phase, fx iigas (T,pi) = + A/i ilgas (T,p;) 

and the free energy of the particle in the adsorbed state 
F i<st = E\°s t - k B T\n(zJ%). To evaluate the desorption 
rate constant using Eq. we thus need to determine, in 
addition to the adsorption rate constant, the vibrational 
partition function of the particle in the adsorbed state, 
zj 1 ^., the temperature and pressure dependent part of the 
gas phase chemical potential, A/ij lgas (T,pj), as well as the 
binding energy, E^ t d = -'X\°t of th e particle at 

its adsorption site st. 



3. 



In a diffusion process a particle i moves from a site st 
to a site st'. If an appropriate saddle point along the 



4- Reaction 

If a suitable reaction coordinate can be identified that 
contains a saddle point, again TST can be used to deter- 
mine the reaction rate constant. In a general form the 
reaction rate constant can then be written as^ 9 - 



' roac 



(T) = 

,rcac.TST m 



(9) 



k B T 



exp 



AE^j 
k B T 



with i denoting the initial state and / the final state 
of the reaction process. As for the diffusion events the 
factor f rea J' T T can be obtained from the ratio of the 
partition functions in the transition and initial state. In 
addition the reaction barrier, AE^j, defined as the en- 
ergy difference in the transition and initial state, has to 
be determined. 



C. Computational Setup 

The most crucial input parameters to determine the 
rate constants as discussed in the previous Section are 
the energy barriers for adsorption, diffusion and reac- 
tion, as well as the binding energies at a certain adsorp- 
tion site. We calculate the total energies that are needed 
to evaluate these quantities using DFT and within the 
full-potential (linearized) augmented plane wave + lo- 
cal orbital (L)APW+lo metho d 31 > 32 as implemented in 
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the WIEN2k code^ 3 -. The exchange-correlation (xc) en- 
ergy is treated within the generalized gradient approxi- 
mation (GGA) using the PBE 34 xc-functional. All sur- 
faces are simulated using the supercell approach with in- 
version symmetric slabs consisting of five Pd(100) layers 
with the reconstructed surface oxide layer plus additional 
O/CO on both sides. We verified that further increas- 
ing the number of Pd(100) layers to seven changes the 
binding energies of both adsorbates by an insignificant 
amount (< 2meV per adsorbate). The vacuum between 
consecutive slabs is at least 14 A. The outermost adsor- 
bate layers, the surface oxide and topmost palladium sub- 
strate layer are fully relaxed. All calculations of gas- 
phase molecules (O, CO, CO2) are done in rectangular 
supercells with side lengths of (13 x 14 x 15) bohr (for O), 
(13xl4xl8)bohr (for 2 and CO) and (13 x 14 x 20) bohr 
(for C0 2 ). 

The muffin-tin radii are set to R^ T = 2.0 bohr for pal- 
ladium, i?MT = 1.0 bohr for oxygen and R^t = 1-0 bohr 
for carbon. Inside the muffin-tins the wave functions 
are expanded up to l^ ax — 12 and the potential up to 
igj^ = 6. For the \/E surface unit cell a [4 x 4 x 1] 
Monkhorst-Pack (MP) grid (8 irreducible k-points) has 
been used to integrate the Brillouin zone (BZ), and for 
larger surface unit cells the MP-grid has been reduced 
accordingly to assure an equivalent sampling of the BZ. 
For the gas-phase molecule calculations in the large rect- 
angular supercells T-point sampling was employed. The 
energy cutoff for the expansion of the wave function in 
the interstitial is E^ l ax = 20 Ry and for the potential 
Emax — 196 Ry. The chosen computational setup is thus 
exactly the same as in the preceding thermodynamic 
study on this system^. As detailed there^ using these 
basis set parameters the reported binding energies per 
oxygen atom or per CO molecule are converged to within 
50meV, which is fully sufficient for the arguments and 
conclusions presented here. 



D. The kMC model 




FIG. 1: (Color online) In the upper panel a schematic top 
view of the v5 surface oxide structure is shown. The large, 
grey spheres represent Pd atoms, the small, red (dark) ones O 
atoms. The Pd atoms of the Pd(100) substrate are darkened. 
The y/E surface unit cell, as well as the (1 x 1) surface unit 
cell of Pd(100) are indicated. As illustrated the y/E surface 
unit-cell consists of two halves that differ only with respect to 
the underlying Pd(100) substrate. We find the difference in 
binding energies at the equivalent adsorption sites in the two 
halves to be very small, and correspondingly use the reduced 
unit-cell in the kMC simulations. In the lower panel these 
prominent adsorption sites and the final kMC lattice are il- 
lustrated. The three high-symmetry adsorption sites (bridge, 
top2f, top4f) are indicated by the yellow (light) crosses. The 
hollow adsorption site of the upper O atoms (red (dark) cross) 
is likewise included as a site in the kMC lattice, whereas the 
reconstructed Pd layer as well as the lower O atoms are fixed. 
In the final kMC lattice only the bridge and hollow sites are 
considered as shown on the lower left (see text). 



Having established the theoretical framework for the 
kMC simulations the next crucial step is to build a suit- 
able lattice model that represents the investigated sys- 
tem. For this, we use the insight gained in our pre- 
ceding constrained thermodynamics study! that identi- 
fied the y/E as most relevant oxidic structure under cat- 
alytically interesting gas phase conditions (pi ~ 1 bar, 
T ~ 300 — 600 K). As initially discussed in Section I wc 
thus concentrate the kMC simulations on this particu- 
lar structure and develop a first-principles kMC lattice 
model that allows for a proper mapping of the y/E sur- 
face oxide structure. 



1. The kMC lattice 

In the upper panel of Fig. Q] a schematic illustration of 
the surface oxide on Pd(100) is shown. The y/E surface 
oxide structure characterized in UHV corresponds essen- 
tially to a (O-Pd-O)-trilayer of PdO(lOl) on top of the 
Pd(100) substrate^ 3 -. Each y/E surface unit-cell contains 
four Pd atoms in the reconstructed oxide layer. In this 
stoichiometric termination of the surface oxide two of the 
Pd atoms are fourfold coordinated and two are twofold 
coordinated by oxygen. There are also two kinds of oxy- 
gen atoms, two of them sit on-top of the reconstructed 
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Pd layer (upper O atoms) and two at the interface to 
the Pd(100) substrate (lower O atoms) forming the pre- 
viously mentioned trilayer structure. 

At first sight the surface oxide structure exhibits six 
high-symmetry adsorption sites within the y5 surface 
unit cell, where either oxygen or CO can be adsorbed,? 
These are two bridge sites (br), two twofold by oxygen 
coordinated top sites (top2f) and two fourfold by oxygen 
coordinated top sites (top4f). The difference between 
each of these doubly existing adsorption sites arises only 
from the arrangement of the underlying Pd(100) lattice. 
These small structural deviations seem to have little con- 
sequences on the binding at these different prominent ad- 
sorption sites. The computed binding energies of O and 
CO in the two bridge sites differ by less than 0.1 eV. Sim- 
ilarly for the two top2f and top4f sites the difference is 
even smaller, Ai? bmd < 50meV (a more detailed discus- 
sion of the binding energetics can be found in Ref. |7|). In 
our kMC model, we therefore neglect these small differ- 
ences and build the lattice from the half sized y5 surface 
unit-cell as shown in Fig. [TJ each containing three dif- 
ferent adsorption sites, bridge, twofold top and fourfold 
top. 

Inspecting the computed energetics at these three dif- 
ferent adsorption sites in more detail the kMC lattice 
model can even be further simplified. We find that in 
the top2f and top4f sites the CO is bound much weaker, 
i.e. by 0.3 — 0.8 eV, compared to the bridge site. In ad- 
dition, the two top sites correspond only to very shallow 
PES minima, with diffusion barriers of less than 50 meV 
to a neighboring bridge site. Considering further that 
due to the computed strongly repulsive interactions^ it is 
not possible to adsorb CO at a top site next to an already 
occupied bridge site, we realize that the metastable top 
sites will effectively never be occupied. For oxygen the 
situation is even more pronounced, since here the addi- 
tional adsorption of an O atom at the y5 surface oxide 
is only possible in the bridge site, whereas the two top 
sites are not even metastable. The two different top sites 
are therefore disregarded in setting up the kMC lattice 
representation of the surface oxide. 

On the other hand, considering only the adsorption 
into prominent high-symmetry sites at the perfect y/5 
structure is not sufficient to address the CO-induced de- 
composition of this lattice. For this, we assume that the 
surface oxide layer starts to destabilize by the removal of 
the upper oxygen atoms from the trilayer structure. We 
therefore also include the corresponding hollow sites that 
are completely occupied by these surface oxygen atoms in 
the perfectly stoichiometric ^/E structure as explicit sites 
in the kMC model, cf. Fig.[Q The basis for the final kMC 
lattice model is thus formed by the fixed reconstructed 
Pd layer and the lower O atom in the PdO trilayer. On 
this lattice there are two different active site types (br 
and hoi, cf. left part of the lower panel in Fig. [T]) . Each 
site can either be vacant, or filled with one oxygen atom, 
or one CO molecule, whereas multiple adsorption at a 
site is not possible. 



Considering all non-concerted processes possible on 
this lattice we obtain a total number of 26 possible pro- 
cesses. This includes five adsorption processes, the uni- 
molecular adsorption of CO into br or hoi sites and the 
dissociative adsorption of O2 into adjacent br-br, br-hol 
or hol-hol sites, as well as the corresponding five des- 
orption processes. Furthermore there are eight site and 
element specific diffusion processes describing the diffu- 
sion of O/CO from br— >br, hoi— > hoi, br— >hol and hoi— >br 
sites, as well as four reaction processes including the re- 
action of O br + CO br , O ho1 + CO ho1 , O br + CO ho1 and 
Qhoi _|_ co br . Since the reaction processes are modeled 
as associative desorption there are an additional four pro- 
cesses for the dissociative adsorption of C02- This list 
of non-concerted processes is exhaustive with respect to 
the chosen kMC lattice setup. However, it can not be 
excluded that there might be also other, more complex, 
concerted processes which could have an influence on the 
here discussed system. 



2. Simulation setup 

All kMC simulations are performed on the final lattice 
model containing (10 x 10) half y/5 surface unit cells, i.e. 
100 br and 100 hoi sites, with periodic boundary condi- 
tions. To ensure that the size of the chosen lattice is suf- 
ficient simulations were also performed on larger lattices 
containing up to 1250 br and hollow sites, respectively, 
without observing significant differences concerning the 
investigated quantities. These quantities of interest for 
the present study are the average surface occupations 
6i tS t of a particle i on a site st under steady-state condi- 
tions, which is evaluated as 

Vi,st = ^ a , > UU) 

where At n is the time step of the nth Monte Carlo move. 
In all simulations we checked carefully using different 
starting configurations and varying simulation lengths 
that the system actually reached steady-state, before the 
average surface occupations were evaluated. 



E. Evaluation of the rate constants 

For all 26 processes that are possible in the kMC 
model, the rate constants have to be determined. The re- 
spective expressions for the four different process types, 
i.e. adsorption, desorption, diffusion and reaction, have 
been described in Section III Bl Here we determine 
the corresponding parameters entering the equations for 
the different rate constants, following the procedure de- 
scribed in more detail in Ref. [29|. 
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1. Adsorption 

To calculate the rate constants of adsorption in Eq. ([3]) 
we specify the active area Ai tSt for the impingement 
of each species on each site as simply half the surface 
area of the surface unit cell in the kMC lattice. As 
shown in Fig. [T] this area corresponds to half the sur- 
face unit cell of the \/5 structure {l/2A^ = 19.47A 2 ) 
i.e. Aco,br = ^4co,hoi = 1/4^4^/5 for the adsorption of 

CO and Ao 2 ,br-br = ^0 2 ,hol-hol = ^0 2 ,br-hol = 1/6A^ 

for the adsorption of O2. 

To obtain an estimate for the local sticking coeffi- 
cient §i s t in Eq. ((3]) we perform DFT calculations to 
determine hrst of all if there are any significant barri- 
ers, AEf^l, along the MEP for any of the five differ- 
ent adsorption processes. For the CO molecule this was 
done by vertically lifting the CO in an upright geome- 
try from a bridge resp. hollow site. For different heights 
of the CO molecule above its adsorption site the energy 
of the whole system was then minimized with respect 
to all other degrees of freedom . For neither of the two 
sites a barrier was found along this desorption resp. ad- 
sorption pathway i.e. AE^ br = AE^ hol = 0.0 eV. 
For the dissociative adsorption of O2 molecules a PES 
was mapped out in which the energy is given as a func- 
tion of the bond length between the two O atoms and 
the height above the surface of their center of mass (so- 
called elbow-plot). The lateral position of the center of 
mass within the surface unit cell, as well as the orienta- 
tion of the O2 molecule were kept fixed. Since there are 
only two hollow sites within the y/E surface unit cell a 
lifting of two oxygen atoms would correspond to a com- 
plete removal of all upper oxygen atoms in the surface 
oxide structure due to the periodic supercell setup of the 
DFT calculations. Therefore, for these calculations the 
size of the \fE surface unit cell was doubled in the direc- 
tion of the involved sites. We find that the dissociation 
of an O2 molecule in this particular orientation appears 
to be non-activated. A similar result is also found for 
the O2 adsorption in neighboring bridge-hollow sites, i.e. 
A£§ d 2 V_ ho] = ^Stbr-ho! = 0.0 eV. For the adsorp- 
tion in two neighboring bridge sites, though, the PES 
reveals a high barrier of A£^ s br _ bl . = 1.9 eV. 

A second parameter entering the sticking coefficient 
in Eq. Q is the factor /f^ 1 representing the reduction 
in sticking due to the fraction of impinging molecules 
not following the MEP. Since the kMC simulations are 
performed in a temperature range of T = 300 — 600 K, 
where the impinging gas phase molecules can still be con- 
sidered as rather slow, it is assumed that the molecules 
are quite efficiently steered along the MEP and thus the 
value of is approximated by /f^f s» 1 for all five 
adsorption processes. With these parameters the stick- 
ing coefficients are unity for the four barrier-free adsorp- 
tion processes. For the activated adsorption of O2 in 
two neighboring bridge sites the sticking coefficient is as 
small as ~ 10~ 16 even for a temperature of T = 600 K. 
Even if the crude way of extracting the approximate acti- 



vation barrier from the restricted elbow plot would thus 
introduce an error of many orders of magnitude in the 
sticking coefficient, this adsorption process will therefore 
still be insignificant compared to the adsorption of O2 in 
two neighboring hollow or neighboring bridge and hollow 
sites. The results of the kMC simulations reported below 
where indeed found to be unaffected, even when lowering 
the activation barrier for this adsorption process by up 
to 0.5 eV. 



2. Desorption 

Having determined the five different adsorption rate 
constants the respective desorption rate constants can 
be calculated using detailed balance via Eq. ([SJ. To 
a first approximation the vibrational partition function 
in the adsorbed state entering Eq. ([5]) is set to unity, 
zj l ^ t w 1, for all five desorption processes. In the 
gas phase this approximation would certainly be rea- 
sonable over the here investigated temperature range 
(T = 300 — 600 K), i.e. only the first vibrational state 
of the O2 and CO molecules is excited. In the adsorbed 
state, though, this might be changed due to the lower 
frequency modes resulting from the adsorbate-substrate 
bond, which could then lead to a somewhat larger par- 
tition function, zj* t > 1, and therefore slightly smaller 
desorption rate constants. Since this enters only linearly 
into the desorption rate constant we nevertheless found 
that considering such increased partition functions did 
not significantly affect the conclusions drawn from the 
simulations below. 

The temperature and pressure dependent part of the 
gas phase chemical potential A^i igas (T,pi) is calculated 
by approximating the O2 and CO gas phase with ideal gas 
phase reservoirs and evaluating the translational, vibra- 
tional, rotational, electronic and nuclear partition func- 
tions using statistical thermodynamics^. This approxi- 
mation by an ideal gas only introduces a negligible error 
in the temperature and pressure range considered tiered. 

A further parameter that is needed to evaluate Eq. ([5]) 
is the binding energy of the particles in their respective 
adsorption sites, Ef 1 ^. . Since the binding energies of O 
and CO in bridge and hollow sites also depend on lateral 
interactions between the adsorbates, we expand them in a 
lattice gas Hamiltonian (LGH)2&i2L22i22i2a. Aiming only 
at a first approximation of the effect of lateral interac- 
tions we employ a crude LGH that is restricted to the 
nearest neighbor pair interactions between adsorbates in 
different and alike sites, yielding 

H = ^[no^+nco,^,,] (11) 

i 

+ / J[Vb-Q,tj n os no,] + Vbo-CO,ij n C o.i n co,j 

ij 

+Vo-CO,ij n O,i n CO,j] , 

where i runs over all lattice sites in the kMC lattice and 
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FIG. 2: (Color online) Considered nearest neighbor pair in- 
teractions between adsorbates on the y/E surface oxide lattice. 
Large, grey spheres represent Pd atoms. Small, red spheres 
oxygen atoms included in the kMC lattice and small yellow 
spheres adsorbed O and/or CO. 



j only over the corresponding nearest neighboring sites. 
Eq i and E^ i are the on-site energies of O and CO 
on a site i, Vij is the interaction parameter between 
two adsorbates sitting on sites i and j and rii is the 
occupation number, i.e. n.; = denotes that site i is 
empty and — 1 that it is occupied. To obtain an es- 
timate of the error introduced by limiting the LGH to 
the nearest neighbor pair interactions only we performed 
DFT calculations with sparser adsorbate arrangements 
in larger surface unit cells, i.e. in (2\/5 x y/5)R27° and 
(VE x 2V5)i?27° cells. The obtained binding energies of 
O and CO in these larger surface unit cells differ by less 
than 0.1 eV from the respective values calculated within 
the y5 surface unit cell. This change in binding energies 
is within the here assumed level of accuracy, rational- 
izing the choice to neglect lateral interactions between 
adsorbates extending beyond the V5 surface unit cell in 
the current study. Employing the LGH in Eq. (fTTj) the 
binding energy of a CO molecule in a site i is then given 
by 

= (12) 
= H(n CO ,i = 1) - H{n CO ,i = 0) 

= EcO.i + 2 ^ [VcO~CO,ij "COj + Vo~CO,ij n O,j] 

j 



TABLE I: Four on-site energies E° and 10 first nearest neigh- 
bor pair interaction parameters V for the lattice gas Hamil- 
tonian describing the adsorption of O and CO in hollow and 
bridge sites on the y/Z surface oxide. All values in eV. 



p0 

^O.br 

-0.51 


-1.95 


^CO.br 
-1.40 


^CO.hol 

-1.92 


Vb-0,br-br 
0.08 


Vb-0,hol-hol 

0.07 


Vb-O.br-hol 

0.08 




VcO-CO,br-br 
0.08 


VcO-CO,hol-hol 

0.13 


VcO-CO.br-hol 
0.14 




Vb-CO,br-br 
0.06 


Vb-CO,hol-hol 

0.11 


Vb-CO,br-hol 
0.13 


Vb-CO,hol-br 

0.12 



\/5 surface unit cell, and expressing the respective energy 
of each configuration in terms of the LGH expansion (cf. 
Eq. (fTTj) ). Since within the kMC model the lower oxygen 

atoms of the V5 surface oxide trilayer form a part of the 
lattice itself the energies of the different configurations 
are calculated with respect to the empty kMC lattice, 
i.e. 



A£(0,CO@v^- 20) 



(13) 



E, 



O,CO'0v / 5-2O 



E 



V5-20 



JV l/2£&? 



NcoEco 



where E^_ 20 corresponds to the total energy of the 
surface oxide structure without the two upper oxygen 



is 



^bind 
O2 



atoms. The total energy of an oxygen molecule, Eq, 
calculated using the relationship Eq* = 2Eq * + 
with a highly converged value of the oxygen gas phase 
binding energy of £^ nd = -6.22 eVii. Table Q] lists the 
determined values for the four on-site energies and 10 in- 
teraction parameters. With these parameters the binding 
energy of every adsorbate for any random configuration 
on the lattice can be evaluated using Eq. (fl2|) . With this 
binding energy the respective desorption rate constants 
can then be calculated via Eq. ([5]). 



Accordingly, the binding energy of an oxygen atom is 
calculated. During the associative desorption of an O2 
molecule, respectively, two neighboring sites are depleted, 
which can be expressed likewise using the LGH. As 
shown in Fig. [2] the first nearest neighbor pair interac- 
tions correspond to interactions between two neighboring 
bridge sites, two neighboring hollow sites and between 
neighboring bridge and hollow sites. This generates 
six different interaction parameters between like species 

(Vb-0,br-br, Vb-O.hol-hol, Vb-0,br-hol, VcO-CO,br-br, 

Vco-co,hoi-hoi, Vco-co,br-hoi) and four different ones 
between unlike species (Vb-co,br-br, Vb-co,hoi-hoi, 
Vb-co,br-hoi, Vb-co,hoi-br)- Together with the four on- 
site energies we thus have to determine 14 parameters for 
the LGH. This is achieved by fitting the parameters to 
DFT total energies of 29 different ordered configurations 
with O and/or CO in bridge and hollow sites within the 



3. 



Within the kMC lattice model for the \fh surface oxide 
there are eight different diffusion processes, four between 
like sites (O/CO hoi— > hoi and br— >br) and four between 
unlike sites (O/CO hoi— >br and br— >hol). To obtain the 
corresponding rate constants of diffusion the vibrational 
partition functions of the adsorbates in the initial state 
and in the transition state are needed, as well as the 
diffusion barriers, cf. Eqs. © and ([7]). 

With a similar justification as for the desorption rate 
constants, we assume the partition functions in initial 
and transition state to be comparable, i.e. z] 



ib 



y vib 



• 1 ,• j.diff,TST 

*«f, ts yielding firtUtf 



1 for all eight diffusion 



processes. To calculate the diffusion barriers, AE' 



diff 
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TABLE II: Barriers for the diffusion between the different 
sites within the \/5 surface unit cell. All values are in eV. 





A P diff 




A rjidiff 


o 

CO 


1.2 
0.4 


1.4 
0.6 


0.1 

0.3 



TS 




final 



initial final initial 

FIG. 3: Schematic illustration of the diffusion barriers be- 
tween energetically like and unlike sites. 

we determine the transition state by defining a reaction 
coordinate connecting the initial and final state and cal- 
culate the energy of O/CO for different positions along 
this coordinate. Since we found the surface oxide trilayer 
to be rather mobile with respect to a registry shift par- 
allel to the Pd(100) substrate^, the lateral positions of 
the reconstructed Pd atoms are additionally kept fixed 
in these calculations. The resulting barriers for these 
approximate diffusion paths are listed in Table [TTJ The 
barriers have been calculated within the ^fh surface unit 
cell with both hollow sites occupied by O (O/CO diffu- 
sion br — > br), only one hollow site occupied by O (O/CO 
diffusion hoi — > br) and one hollow site occupied by ei- 
ther O or CO (O/CO diffusion hoi -► hoi). The initial 
and the final states of a diffusion process can either be 
energetically degenerate (left part of Fig. [3]) or - for the 
diffusion between unlike sites or like sites with a differ- 
ent nearest neighbor coordination - the energies of the 
initial and final state wells differ (right part of Fig. [3} . 
In the latter case the values listed in Table |TT] represent 
the minimum value of the diffusion barrier, i.e. if the 
final state has the same or a lower energy than the ini- 
tial state the diffusion barrier is equal to the tabulated 
value. If on the other hand the final state is higher in 
energy the difference in binding energies between the ini- 
tial and final state is added to the listed barrier to ensure 
that the detailed balance criterion is fulfilled for the two 
time-reversed diffusion processes. Since this difference in 
binding energies depends on the nearest neighbor coor- 
dination of the two lattice sites involved in the diffusion 
process, lateral interaction effects are this way also con- 
sidered in the diffusion process rate constants. 

The described approximations in determining the ex- 
act transition states lead to uncertainties in the stated 
barrier values which in turn influence the corresponding 
diffusion rate constants. We checked on the importance 
of the various diffusion processes on the simulation results 
and found them to play only a minor role. Essentially, 



TABLE III: Barriers for the reaction of O and CO adsorbed 
in different sites within the v5 surface unit cell. All values 
are in eV. 



O br +CO br 


Qh0l +C Qh0l 


br + C Qhol Qhol +C Qbr 


A£ reac 1.0 


1.6 


0.5 0.9 



even completely switching off all diffusion processes led 
only to insignificant changes of the results presented be- 
low. We are thus confident that the rather coarse proce- 
dure to determine the diffusion barriers is sufficient with 
respect to the current purpose of the simulations. 



4- Reaction 

To calculate the rate constants for the four possible 
reaction processes (O br +CO br , O hol +CO ho \ O br +CO ho1 
and O hol +CO br ) the factor f™ c / TST w 1 in Eq. © is 
approximated by unity. Again it is thus assumed that to 
a first approximation the partition functions of the initial 
and transition state are in the same order of magnitude, 
and the ratio of the two is thus close to one. 

To determine the reaction barriers, AE?fJ5, the cor- 
responding transition states are located by dragging the 
two reactants towards each other and mapping out the 
corresponding 2D-PES (O br +CO ho1 and O hol +CO br ). 
Hereby only the coordinates of the reactants along the 
dragging direction are kept fixed whereas all other coor- 
dinates are fully relaxed. Similar to the problem encoun- 
tered when determining the diffusion barriers the recon- 
structed palladium layer showed a rather high mobility 
with respect to a lateral shift along the substrate follow- 
ing the dragging direction. To circumvent this problem 
a two-step approach was chosen to nevertheless identify 
appropriate transition states. In the first step the lateral 
coordinates of the reconstructed palladium layer were ad- 
ditionally fixed and the energy is only minimized with re- 
spect to the remaining coordinates. In a second step the 
identified transition state geometry is relaxed with re- 
spect to the lateral coordinates of the topmost Pd layer 
while fixing the coordinates of the reactants O and CO, 
to release some of the lateral stress that was induced by 
fixing the lateral positions of the Pd atoms in the first 
step. For the reaction of O br +CO br and O hol +CO ho1 this 
two-step approach was not necessary. Similar to the dis- 
sociative adsorption of O2 an elbow plot was calculated 
giving the energy as a function of the O-CO distance 
and the height of the corresponding center of mass above 
the surface. Due to the symmetry of the calculated ge- 
ometries no shift of the substrate layer could occur in this 
case. All correspondingly calculated reaction barriers are 
summarized in Table [TTT1 

For the reaction processes the energies of the initial, 
transition and final states are also influenced by the near- 
est neighbor interactions. Inspecting the geometries of 
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the transition states they appear to be more related to 
the initial states and we therefore assume that the tran- 
sition states are similarly affected by these interactions 
as the initial states. Thus, the energy difference between 
the initial and transition states is not influenced by the 
nearest neighbor interactions and the reaction barriers 
are constant with respect to different configurations on 
the surface. 

Since the reaction processes are modeled as associative 
desorption the corresponding time reversed process is the 
dissociative adsorption of CO2 into the respective sites. 
The adsorption barrier for such a process can be obtained 
from the reaction barrier and the corresponding binding 
energies of O and CO on the surface with respect to the 
free CO2 molecule. The lowest of the four resulting bar- 
riers for the dissociative adsorption of CO2 (leading to 
CO in bridge and O in a neighboring hollow site) is still 
very large with AE^q 2 br _ ho j = 1.4 eV. At a temperature 

of T = 600 K the local sticking coefficient 5co 2 ,br- hoi is 
even in this case then only of the order of 10 -12 . Even 
if the CO2 generated at the surface is not readily trans- 
ported away and a noticeable CO2 pressure would build 
up close to the surface, the re-adsorption of CO2 on the 
surface is therefore still negligible due to the very low 
sticking coefficient and is thus not explicitly considered 
in the kMC simulations. 



III. RESULTS 

Having established the setup of the kMC lattice model 
and the respective rate constants for the 26 different pro- 
cesses based on DFT energetics and harmonic TST we 
proceed to the results of the kMC simulations. Although 
the established kMC model has a wider range of applica- 
bility, we focus here on the stability of the thin surface 
oxide layer during steady-state catalytic CO oxidation at 
the surface. 



A. Reproducing the thermodynamic "phase" 
diagram 

In order to make contact to our preceding constrained 
thermodynamics worki, we perform the kMC simulations 
in a first step only considering adsorption, desorption and 
diffusion processes, whereas the reaction processes are ex- 
cluded. The occupation of the different sites at the sur- 
face under steady-state conditions reflects then the same 
situation as a surface in a constrained thermodynamic 
equilibrium with an oxygen and CO gas-phase, but in- 
cluding the effect of configurational entropy 6 . The results 
of the constrained thermodynamics approach arc summa- 
rized in Fig. 2] (for a detailed discussion see Ref.0). In 
the shown "phase" diagram the most stable surface struc- 
tures for any given chemical potential of the O2 and CO 
gas-phase are presented (the respective pressure scales for 
T = 300 K and 600 K are given in the upper x-axes and 
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FIG. 4: (Color online) Thermodynamic surface "phase" dia- 
gram of the Pd(100) surface in a constrained equilibrium with 
and O2 and CO gas-phase as described in detail in Ref. 0. 
The differently shaded areas mark the stability regions of var- 
ious surface structures for a given chemical potential of the O2 
and CO gas-phase. Blue/solid areas include O/CO adsorp- 
tion structures on Pd(100), red/ white-hatched areas comprise 
surface oxide structures and the grey/crosshatched area rep- 
resents the stability range of PdO bulk oxide. For convenience 
the dependence of the chemical potential of the two gas-phases 
is translated into pressure scales for T = 300 K and T = 600 K 
(upper x axes and right y axes). For the stability region of the 
surface oxide structures additionally schematic figures of the 
different surface configurations are shown. Grey, large spheres 
represent Pd atoms, red (dark) small spheres oxygen atoms 
and yellow (light) small spheres CO molecules. The white ar- 
rows indicate the pressure conditions of the kMC simulations 
for a temperature of T = 300 K, T = 400 K and T — 600 K. 
The oxygen pressure is always fixed at po 2 — 1 atm, while the 
CO pressure is varied between I0~ < pco < 10 atm. 



right y-axes). The stability regions of the different sur- 
face structures are denoted by differently shaded areas. 
For the here discussed kMC simulations it is particularly 
important that the different "phases" can be subdivided 
into three groups. All areas colored in blue/solid rep- 
resent structures with O and CO adsorbed on the pris- 
tine Pd(100) surface, all red/ white-hatched areas com- 
prise structures involving the reconstructed y/E surface 
oxide, and the grey/crosshatched area marks the stabil- 
ity region of the bulk oxide. As discussed in Section Hi Dl 
the here chosen kMC lattice model is based on the lat- 
tice structure of the y/b surface oxide and does there- 
fore not allow for changes between this structure and the 
Pd(I00) surface or the bulk oxide. Consequently, the con- 
strained kMC simulations can only reproduce the part of 
the "phase" diagram where structures based on the sur- 
face oxide are most stable, i.e. the red/white-hatched 
region in Fig. [U 

Within the constrained kMC simulations we now eval- 
uate the steady-state average occupation of bridge and 
hollow sites and compare the obtained structures to the 
structures predicted within the thermodynamic surface 
"phase" diagram. As mentioned above the two ap- 
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proaches should yield equivalent results, except at the 
"phase" boundaries where configurational entropy is ex- 
pected to become important. The simulations are per- 
formed for (T, p)-conditions that cover the entire stability 
region of the surface oxide (red/white- hatched region in 
Fig. HI). Specifically, this means for T — 300 K pressures 
of po 2 = W- w - 10 10 atm and p C o = 10~ 10 - latm, 
and for T — 600 K pressures of po 2 = 1 — 10 10 atm and 
Pco = 1 — 10 6 atm, respectively. We find that the kMC 
simulations reproduce accurately all four different surface 
oxide "phases" that appear in this part of the constrained 
"phase" diagram. In all four "phases" all hollow sites are 
always filled with oxygen. The bridge sites are either 
empty, half filled with O or CO, or completely filled with 
CO depending on the respective gas-phase conditions and 
as detailed by the schematic representations in Fig. 2] 
As expected from the common underlying DFT energet- 
ics, we therefore find the results of the constrained kMC 
simulations to be fully consistent with those obtained 
previously within the constrained thermodynamics ap- 
proach. The more extended sampling of configurations 
due to the LGH expansion within the kMC simulations 
does not lead to the appearance of new, stable ordered 
structures in the constrained "phase" diagram and thus 
confirms that the relevant configurations were included 
in the preceding thermodynamic approach. The effect of 
configurational entropy in the constrained kMC simula- 
tions is thus verified to "only" smear out the transitions 
between the different ordered "phases" at the finite tem- 
peratures considered^, but does not affect the location of 
the "phase" boundaries themselves. Due to the hereby 
confirmed consistency between the two approaches, we 
can directly compare the results in the next Section, when 
allowing the reaction events to occur in the kMC simu- 
lation. This brings us in the position to directly address 
the effect of the kinetics of the ongoing chemical reactions 
on the surface structure and composition. 



B. Onset of surface oxide decomposition under 
reaction conditions 

By allowing the four different reaction processes to oc- 
cur in the kMC simulations, we now proceed to evaluate 
the stability of the surface oxide structure under reac- 
tion conditions. Since the kMC lattice is fixed to the 
structure of the y/5 surface oxide, the explicit structural 
change corresponding to the transition from the surface 
oxide to a pristine Pd(100) surface as predicted in the 
thermodynamic "phase" diagram (red/white- hatched ar- 
eas — > blue solid areas in Fig. [4]) cannot directly be ad- 
dressed. Instead, we start our simulations always in a 
(T, p)-rangc where the surface oxide is definitely the most 
stable "phase" (high oxygen and low CO chemical po- 
tential), and where our kMC lattice model is certainly 
valid. Keeping the oxygen pressure and the temperature 
fixed we then increase the CO pressure as indicated by 
the white, vertical arrows in Fig. 2] Exceeding a certain 



pressure ratio of Po 2 /Pco, we expect the onset of the de- 
composition of the y5 structure to be characterized by 
a depletion of the otherwise always essentially fully with 
oxygen occupied hollow sites. In the kMC simulations 
we therefore monitor the average occupation of hollow 
sites by oxygen atoms, 6>o tol > during steady-state. If all 
hollow sites are occupied by oxygen, i.e. #o ho1 ~ 100%, 
then the overall surface composition clearly represents 
the intact surface oxide structure (this is the case in all 
four "phases" of the surface oxide structure in Fig. 0|. 
If on the other hand at increasing pco oxygen atoms 
get increasingly depleted beyond the amount of thermal 
fluctuations that create isolated, uncorrelated vacancies 
or double vacancies due to oxidation reaction or oxygen 
desorption events, then the deficiency of O atoms at the 
surface might destabilize the v5 surface structure and 
lead to a local lifting of the surface oxide reconstruc- 
tion. Here, we therefore interpret a drop of the average 
occupation below 90% as the onset of a destabilization 
of the surface oxide. Since the average occupation is a 
global measure over the whole simulation cell this quan- 
tity would be a bad indicator, if there was an appreciable 
effective attraction between formed O vacancies. A de- 
crease in the average occupation of hollow sites by oxygen 
of 10% could then involve a complete, local depletion of 
10% of the surface oxide area. In our kMC simulations 
such a locally concentrated depletion of oxygen atoms is 
not observed. Instead, we find the formed oxygen vacan- 
cies to be distributed over the entire simulation cell and 
are thus confident that monitoring the average occupa- 
tion, O hoi, is a suitable measure to determine the onset 
of surface oxide decomposition. 

The simulations are performed for three different tem- 
peratures, T = 300,400 and 600 K. For each tempera- 
ture the oxygen pressure is fixed to po 2 — 1 a t m & n d the 
simulations are run for different CO pressures between 
pco = 10~ 5 — 10 5 atm covering the gas-phase conditions 
marked by the white arrows in Fig. 2J In Fig. [5] the 
corresponding simulation results are summarized. The 
steady-state average occupation of hollow sites with oxy- 
gen, 0Qhai, is plotted vs. the CO pressure. The solid 
(green) lines represent the results obtained by using the 
calculated reaction barriers as listed in Table IIIII As 
has been discussed in Section UH several, partly crude 
approximations were employed in the determination of 
the individual rate constants entering the kMC simula- 
tions. Before we proceed to analyze the obtained simu- 
lation results, we therefore critically assess the propaga- 
tion of these underlying approximations and uncertain- 
ties to the final kMC results. We checked on this error 
propagation by systematically varying the rate constants 
over several orders of magnitude (corresponding to varia- 
tions in the respective energetic barriers of ±0.2 eV) and 
each time monitoring the effect on the average occupa- 
tion displayed in Fig. [5] In almost all cases, variations 
of the rate constants had an insignificant effect, which 
as already described concerns in particular all diffusion 
processes. Only if altering the barrier for the reaction 
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FIG. 5: (Color online) Average occupation of hollow sites 
by oxygen vs. CO pressure for T = 300 K, T = 400 K and 
T = 600 K. The vertical black line marks the boundary be- 
tween the surface oxide and a CO covered Pd(100) surface 
as determined within the constrained thermodynamics ap- 
proach. The influence of the barrier for the O hol +CO br reac- 
tion is illustrated by showing the kMC simulation results for 
two barrier values, 0.9 eV and 0.8 eV (see text). 



process of an oxygen in a hollow site with a CO in a 
neighboring bridge site a noticeable change in the aver- 
age surface populations occurs. To get an estimate for 
the corresponding uncertainties in our presented results, 
we therefore include in Fig. O additionally the results for 
simulations using a 0.1 eV lower barrier for this reaction 
process O hol +CO br , i.e. A££ h a c 1+C0br = 0.8 eV, while all 
other barriers have been left unchanged (dashed (orange) 
lines) . The differences between these results (dashed (or- 



ange) line) and the results for the unmodified barrier 
(solid (green) line) provide then an indication of the un- 
certainty underlying our approach. 

Discussing first the top graph in Fig. [5] it can be seen 
that for a temperature of T = 300 K and an oxygen 
pressure of po 2 = 1 atm the surface oxide is clearly sta- 
bilized for CO pressures up to pco < 10 atm, i.e. 
under such gas-phase conditions 9 hai > 95%. If the 
CO pressure is further increased the O population at 
hollow sites starts to decrease and for CO pressures of 
Pco > 1 atm the surface oxide is certainly completely 
destabilized (0 O h °' ~ 0%). For the slightly lower reac- 
tion barrier (dashed (orange) line) a steeper decrease in 
the oxygen population is found, and the surface oxide is 
stabilized for CO pressures of pco < 10~ 2 atm. To com- 
pare these results to the ones obtained within the con- 
strained thermodynamics approach the stability region of 
the surface oxide as determined in the "phase" diagram 
in Fig. [4] is indicated by the vertical black lines in the 
three graphs of Fig. [5] Following the white arrow in the 
"phase" diagram in Fig.|3]for T — 300 K and po 2 = 1 atm 
this boundary between the surface oxide structure and a 
CO covered Pd(100) surface is reached for a CO pres- 
sure of pco ~ 2 atm. Within the constrained thermody- 
namics approach the surface oxide will thus be stable for 
Pco ^ 2 atm, whereas for pco > 2 atm the CO covered 
Pd(100) surface is the most stable "phase". Comparing 
this to the kMC simulations which include the effect of 
the CO2 formation kinetics, we see that the depletion 
of oxygen atoms in hollow sites and thus the destabi- 
lization of the surface oxide structure starts already at 
slightly lower CO pressures compared to the transition 
from the stability region of the surface oxide to the one 
of a CO covered Pd(100) surface in the "phase" diagram. 
Although we find the kMC simulations to be dominated 
by adsorption/desorption processes of the reactants (as 
anticipated in the constrained thermodynamic approach 
assuming equilibrium with the gas-phase reservoirs), the 
kinetics of the ongoing oxidation reaction thus still leads 
to a slight decrease in the stability of the surface oxide. In 
turn the kMC simulations predict that at a temperature 
of T = 300 K oxygen-rich conditions (pco/po 2 ~ 1 : 10) 
are needed to stabilize the surface oxide structure. 

In Fig. Elb) corresponding results are shown for T = 
400 K. Compared to the simulations at T — 300 K much 
more reaction processes take place in the kMC simula- 
tions at this higher temperature. For CO pressures of 
Pco > 10 atm though almost no reaction processes are 
detected in the steady-state, since under these conditions 
the surface is already completely filled with CO. Then 
again this corresponds to gas-phase conditions where 
also in the constrained thermodynamics approach a CO 
covered Pd(100) surface is the most stable structure 
(cf. Fig. ij). Similar to the results for T = 300 K the 
surface oxide structure is completely stabilized for CO 
pressures of pco < 10 _1 atm. At a pressure ratio of 
Pco/po 2 —2:1, i.e. here pco = 2 atm, there are 
still 94% of all hollow sites occupied by oxygen (solid 
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(green) line in Fig.[5lb)). This result depends, however, 
critically on the barrier for the reaction of O hol +CO br . 
If this barrier is decreased by only 0.1 eV the occupa- 
tion drops to 70%, and at an even lower barrier of 
Ai?Qho C i +CO br = 0.7 eV 0Q ho1 decreases even further to 
~ 30 % under these conditions. Whether or not the sur- 
face oxide prevails at this temperature for the catalyti- 
cally most relevant stoichiometric partial pressure ratio 
can thus not be safely concluded in light of the uncer- 
tainties of our approach. For CO pressures lower than 
Pco < 10 _1 atm, however, the surface oxide is definitely 
stabilized, just as much as it is definitely destabilized at 
CO pressures of pco > 10 atm. We can therefore deter- 
mine the onset of the surface oxide decomposition to fall 
into the pressure range 0.1 atm < pco < 10 atm. Even 
considering the uncertainties in the reaction rate con- 
stants and the limitation to a single kMC lattice, these 
results agree thus rather nicely with the reactor STM 
experiments by Hendriksen et al^ which have been per- 
formed under similar temperature and pressure condi- 
tions (T = 408K, ptot = Pco+Po 2 = l-23atm). Depend- 
ing on the pressure ratio of oxygen and CO Hendriksen 
et al. observed a change in the morphology of the sur- 
face, which was assigned to a change from an oxidic to 
a reduced surface structure and oxygen-rich feeds were 
required to stabilize the oxidic structure. 

The simulation results for the highest considered tem- 
perature of T = 600 K are shown in Fig. [51(c). We find 
that for this elevated temperature the surface oxide is 
now actually stable up a rather sizable CO pressures. 
Only for CO pressures as high as pco = 10 atm the sur- 
face oxide starts to decompose (0 O hoi < 95%), which al- 
ready coincides almost with the border of its stability as 
determined in the constrained thermodynamics approach 
(black vertical line) . In the "phase" diagram in Fig. 2] it 
can be seen that for a temperature of T = 600 K and an 
oxygen pressure of po 2 = 1 atm CO is not stabilized in 
large amounts as adsorbate on the surface oxide. The 
CO is thus not readily available as adsorbed reaction 
partner and the stability of the surface oxide is hardly 
affected by the ongoing catalytic CO2 formation. Up 
to the stoichiometric pressure ratio of pco /po 2 =2:1 
no decomposition of the surface oxide structure occurs, 
and as apparent from Fig. \5j[c) this result holds even in 
light of the estimated uncertainty of the underlying DFT 
energetics. Comparing the critical pco/po 2 ratio deter- 
mined for the decomposition onset in the temperature 
range T = 300 — 600 K, we can thus clearly identify an 
increasing stability of the surface oxide with increasing 
temperature, which at the highest temperatures studied 
reaches well up to the catalytically most relevant feeds. 
Furthermore we find that under these feeds the simulated 
turnover frequencies for the intact surface oxide alone are 
already of similar order of magnitude as those reported 
by Szanyi and Goodman^ for the Pd(100) surface under 
comparable gas-phase conditions. While a quantitative 
comparison is outside the scope of the present work, we 
note that contrary to the prevalent general preconception 



at least this particular surface oxide is thus clearly not 
"inactive" with respect to the oxidation of CO. 



IV. SUMMARY 

We have employed first-principles kMC simulations to 
investigate the stability of the y/E surface oxide struc- 
ture at Pd(100) against CO-induced decomposition un- 
der steady-state operation conditions of the CO oxida- 
tion reaction. The focus has been particularly set on 
the thin surface oxide, since this structure was identi- 
fied as the most stable oxidic "phase" under catalytically 
relevant gas-phase conditions in a preceding constrained 
first-principles thermodynamics study-i. The developed 
kMC model describes the CO oxidation on the \f% surface 
oxide, accounting for all non-concerted processes that are 
possible on the chosen kMC lattice. This includes 26 site 
and element specific adsorption, desorption, diffusion and 
reaction processes. The respective kMC rate constants 
have been obtained from DFT total energy calculations 
together with harmonic transition state theory. 

Monitoring the onset of the surface oxide decompo- 
sition at increasing CO pressures, the central result of 
the kMC simulations is that much in contrast to thicker 
bulk-like oxide films the stability of the V5 surface oxide 
can well extend to gas-phase conditions that are relevant 
for technological CO oxidation catalysis. While the ac- 
curacy of present-day DFT functionals does not permit 
us to conclude on precise values for the critical CO:02 
partial pressure ratio above which decomposition sets in, 
our results show undoubtedly that at ambient pressures 
this ratio is close to the most interesting stoichiometric 
feed conditions (pco : Po 2 =2:1) and will furthermore 
shift decisively towards more CO-rich conditions with in- 
creasing temperature. 

On the basis of these results we conclude that the sur- 
face oxide structure can actually be present under cat- 
alytic reaction conditions at ambient pressures and our 
simulated turnover frequencies indicate that it is then 
also catalytically active towards the oxidation of CO. 
For the catalytically most relevant feeds the system ap- 
pears furthermore to be rather close to a transition be- 
tween the yf& surface oxide structure and a reactant cov- 
ered Pd(100) surface. Particularly in stoichiometric to 
slightly O-rich feeds even local fluctuations in both tem- 
perature or the reactant partial pressures might then lead 
to significant local or global changes in the structure and 
composition of the surface. Under steady-state operation 
conditions oscillations between these two states are thus 
conceivable, so that in its active state the catalyst surface 
might exhibit parts that correspond to Pd(100) covered 
by the reactants, as well as by patches of surface oxide, 
and potentially it is precisely the on-going formation and 
decomposition of the oxidic structure that is key to un- 
derstand the catalytic function of this surface under such 
technologically relevant pressures. For a full understand- 
ing of the catalytic CO oxidation at the Pd(100) surface 
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in corresponding environments it thus appears not to be DFG for support within the priority program SPP1091. 
sufficient to concentrate on either the reduced Pd(100) 
surface or on an oxidic structure, but both surface states 
as well as the morphological transition between the two 
will have to be included in an appropriate model. 
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